
## load data
surv <- read.csv("data/lucid experimental data.csv", colClasses = "character")[c(-1,-2),]

# some respondents did not take the survey after reading consent
surv <- surv[surv$Q1 == "Yes",]

# some respondents never finished
# table(surv$T1Q1 == "" & surv$T2Q1 == "")
# table(surv$T1Q2 == "" & surv$T2Q2 == "")
drops <- surv$T1Q1 == "" & surv$T2Q1 == "" & surv$T1Q2 == "" & surv$T2Q2 == ""
surv <- surv[drops == FALSE,]

## rename questions
names(surv)[c(
   2:7,
   8,9,
   14:16,
   17:19)] <- 
  c("thrm_cong","thrm_corp","thrm_ent","thrm_cath","thrm_who","thrm_nfl",
   "party","ideology",
   "med_job","med_inq","med_pol",
   "employment","emp_size","emp_size_unemp")
head(surv)

## reassign answers for repeated questions to final answer
# treatment indicator variable
surv$treatment <- "1"
surv$treatment[surv$T2Q1 != ""] <- "2"

# support for trade and trade agreements
surv$TradeOpenFavor <- NA
surv$TradeOpenFavor[surv$treatment == "1"] <- surv$T1Q1[surv$treatment == "1"]
surv$TradeOpenFavor[surv$treatment == "2"] <- surv$T2Q1[surv$treatment == "2"]

surv$TradeDealFavor <- NA
surv$TradeDealFavor[surv$treatment == "1"] <- surv$T1Q2[surv$treatment == "1"]
surv$TradeDealFavor[surv$treatment == "2"] <- surv$T2Q2[surv$treatment == "2"]

# employer size
surv$emp_size_comb <- surv$emp_size; surv$emp_size_comb[surv$emp_size_comb == ""] <- surv$emp_size_unemp[surv$emp_size_comb == ""]

# race
race2 <- rep("Other non-white", nrow(surv))
race2[surv$ethnicity == 1] <- "White"
race2[surv$ethnicity == 2] <- "Black"
race2[surv$hispanic %in% c(1,15) == FALSE] <- "Latino"
race2[surv$hispanic %in% c(3:14)] <- "AAPI"
surv$race <- factor(race2, levels = c("White","Black","Latino","AAPI","Other non-white"))

# gender
gender2 <- rep("Female", nrow(surv))
gender2[surv$gender == 1] <- "Male"
surv$gender <- factor(gender2, levels = c("Female","Male"))

## create numeric versions of variables
# support for trade and trade agreements
surv$TradeOpenFavorNum <- NA; surv$TradeOpenFavorNum[surv$TradeOpenFavor == "Oppose a great deal"] <- 1
surv$TradeOpenFavorNum[surv$TradeOpenFavor == "Oppose somewhat"] <- 2; surv$TradeOpenFavorNum[surv$TradeOpenFavor == "Neither favor nor oppose"] <- 3
surv$TradeOpenFavorNum[surv$TradeOpenFavor == "Favor somewhat"] <- 4; surv$TradeOpenFavorNum[surv$TradeOpenFavor == "Favor a great deal"] <- 5

surv$TradeDealFavorNum <- NA; surv$TradeDealFavorNum[surv$TradeDealFavor == "Oppose a great deal"] <- 1
surv$TradeDealFavorNum[surv$TradeDealFavor == "Oppose somewhat"] <- 2; surv$TradeDealFavorNum[surv$TradeDealFavor == "Neither favor nor oppose"] <- 3
surv$TradeDealFavorNum[surv$TradeDealFavor == "Favor somewhat"] <- 4; surv$TradeDealFavorNum[surv$TradeDealFavor == "Favor a great deal"] <- 5

# mediators convert to numeric.
surv$med_job <- factor(surv$med_job, levels = c("Strongly disagree","Somewhat disagree","Neither agree nor disagree","Somewhat agree","Strongly agree"))
surv$med_job <- as.numeric(surv$med_job)
surv$med_inq <- factor(surv$med_inq, levels = c("Strongly disagree","Somewhat disagree","Neither agree nor disagree","Somewhat agree","Strongly agree"))
surv$med_inq <- as.numeric(surv$med_inq)
surv$med_pol <- factor(surv$med_pol, levels = c("Strongly disagree","Somewhat disagree","Neither agree nor disagree","Somewhat agree","Strongly agree"))
surv$med_pol <- as.numeric(surv$med_pol)

# party and ideology, turn to numeric scales
surv$party <- factor(surv$party, levels = c("Strong Democrat","Weak Democrat","Independent - leaning Democrat",
  "Independent - no preference","Independent - leaning Republican","Weak Republican","Strong Republican"))
surv$party <- as.numeric(surv$party)
surv$ideology <- factor(surv$ideology, levels = c("Extremely liberal","Liberal","Slightly liberal","Moderate; middle of the road",
  "Slightly conservative","Conservative","Extremely conservative"))
surv$ideology <- as.numeric(surv$ideology)
surv$income <- as.numeric(surv$hhi)
surv$income[surv$income == -3105] <- NA

# age
surv$age <- as.numeric(surv$age)

# collapse education to numeric binary
surv$college <- as.numeric(surv$education %in% c(4,5,6,7,8))
surv$college[surv$college == -3105] <- NA

# collapse employed into smaller number of categories and make numeric
surv$employed <- as.numeric(surv$employment %in% c("Employed full time","Employed part time"))
surv$unemployed <- as.numeric(surv$employment %in% c("Unemployed looking for work","Unemployed not looking for work"))
surv$retstudis <- as.numeric(surv$employment %in% c("Retired","Student","Disasbled"))

surv$emp_size_comb <- factor(surv$emp_size_comb, c("1-5","6-19","20-49","50-199","200-999","1,000-9,999","More than 10,000"))
surv$emp_size_comb <- as.numeric(surv$emp_size_comb)

# convert empsize into a above/below median dummy
surv$emp_size_big <- as.numeric(as.numeric(surv$emp_size_comb) > median(as.numeric(surv$emp_size_comb), na.rm = TRUE))

# convert used thermometer into above/below median dummies.
surv$thrm_corp_pos <- as.numeric(as.numeric(surv$thrm_corp) > median(as.numeric(surv$thrm_corp), na.rm = TRUE))

## create cleaned version of the data for analysis
surv <- surv[,c("TradeOpenFavorNum","TradeDealFavorNum",
  "treatment",
  "age","gender","race",
  "college","income","employed","unemployed",
  "party","ideology",
  "thrm_corp","thrm_corp_pos","emp_size_comb","emp_size_big",
  "med_job","med_inq","med_pol"
)]

num_vars <- c("TradeOpenFavorNum","TradeDealFavorNum",
  "age",  
  "college","income","employed","unemployed",
  "party","ideology",
  "thrm_corp","thrm_corp_pos","emp_size_comb","emp_size_big",
  "med_job","med_inq","med_pol")
surv[,num_vars] <- apply(surv[,num_vars], 2, as.numeric)

##############################
## Average treatment affect ##
##############################

mod1 <- lm(TradeOpenFavorNum ~ I(treatment=="1"), data = surv); 
mod2 <- lm(TradeOpenFavorNum ~ I(treatment=="1") + age + gender + race, data = surv); 
mod3 <- lm(TradeOpenFavorNum ~ I(treatment=="1") + age + gender + race + college + income + employed + unemployed, data = surv); 
mod4 <- lm(TradeOpenFavorNum ~ I(treatment=="1") + age + gender + race + college + income + employed + unemployed + 
  party + ideology, data = surv); 

mods <- list(summary(mod1)$coef, summary(mod2)$coef, summary(mod3)$coef, summary(mod4)$coef)
vars <- c(rownames(mods[[4]])[2:14], "(Intercept)")
tab <- matrix(data = NA, ncol = length(mods), nrow = length(vars)*2)
rownames(tab) <- c(rep(vars, each = 2))
for(i in 1:length(mods)){mod <- mods[[i]]; tabinmod <- unique(rownames(tab)[rownames(tab) %in% rownames(mod)])
  stars <- rep("", length(tabinmod)); stars[mod[tabinmod,4] < .1] <- "^{+}"; stars[mod[tabinmod,4] < .05] <- "^{*}"; 
  stars[mod[tabinmod,4] < .01] <- "^{**}"; stars[mod[tabinmod,4] < .001] <- "^{***}";  
  tab[rownames(tab) %in% rownames(mod),i] <- c(rbind(paste(myround(mod[tabinmod,1],2), stars, sep = ""), paste("(", myround(mod[tabinmod,2],2), ")", sep = "")))}
rownames(tab) <- c(rbind(vars, "")); rownames(tab)[c(1,3,5,7,9,11,13,15,17,19,21,23,25,27)] <- 
  c("Treated (Large firms ben./Small firms harmed)","Age","Male","Black","Latino","AAPI","Other non-white","College-educated","Income","Employed","Unemployed",
    "Party (D=1,R=7)","Ideology (L=1,C=7)","Intercept") 
N <- c(sum(summary(mod1)$df[1:2]),sum(summary(mod2)$df[1:2]),sum(summary(mod3)$df[1:2]),sum(summary(mod4)$df[1:2]))
N <- paste("\\multicolumn{1}{c}{", N, "}", sep = ""); tab <- rbind(tab, N)
addtorow <- list(); addtorow$pos <- list(28,29); addtorow$command <- c(paste0('\\midrule '),paste0('\\midrule '))
ptable(cbind(rownames(tab), tab), "paper/tables/TableB2.tex")

# CI's
LBs <- c(coef(mod1)[2] - qnorm(.975)*summary(mod1)$coef[2,2],
  coef(mod2)[2] - qnorm(.975)*summary(mod2)$coef[2,2],
  coef(mod3)[2] - qnorm(.975)*summary(mod3)$coef[2,2],
  coef(mod4)[2] - qnorm(.975)*summary(mod4)$coef[2,2])
UBs <- c(coef(mod1)[2] + qnorm(.975)*summary(mod1)$coef[2,2],
  coef(mod2)[2] + qnorm(.975)*summary(mod2)$coef[2,2],
  coef(mod3)[2] + qnorm(.975)*summary(mod3)$coef[2,2],
  coef(mod4)[2] + qnorm(.975)*summary(mod4)$coef[2,2])
CIs <- apply(cbind(LBs,UBs), 2, myround, 2)

est <- tab[1,]
ci <- paste("\\multicolumn{1}{c}{", paste("$[", CIs[,1], ",", CIs[,2], "]$\\phantom{'''''}", sep = ""), "}", sep = "")
ci <- gsub("-","\\\\sminus", ci); ci <- gsub("0\\.","\\.", ci);
N <- c(sum(summary(mod1)$df[1:2]),sum(summary(mod2)$df[1:2]),sum(summary(mod3)$df[1:2]),sum(summary(mod4)$df[1:2]))
N <- paste("\\multicolumn{1}{c}{", N, "\\phantom{''''}}", sep = ""); tab <- rbind(tab, N)
tab <- rbind(est, ci, N)
rownames(tab) <- c("Average treatment effect", "ATE 95\\% CI", "N")
addtorow <- list(); addtorow$pos <- list(1); addtorow$command <- c("")
ptable(cbind(rownames(tab), tab), "paper/tables/TableB1_top.tex")

# 
mod1 <- lm(TradeDealFavorNum ~ I(treatment=="1"), data = surv); summary(mod1)
mod2 <- lm(TradeDealFavorNum ~ I(treatment=="1") + age + gender + race, data = surv); summary(mod2)
mod3 <- lm(TradeDealFavorNum ~ I(treatment=="1") + age + gender + race + college + income + employed + unemployed, data = surv); summary(mod3)
mod4 <- lm(TradeDealFavorNum ~ I(treatment=="1") + age + gender + race + college + income + employed + unemployed + 
  party + ideology, data = surv); summary(mod4)

mods <- list(summary(mod1)$coef, summary(mod2)$coef, summary(mod3)$coef, summary(mod4)$coef)
vars <- c(rownames(mods[[4]])[2:14], "(Intercept)")
tab <- matrix(data = NA, ncol = length(mods), nrow = length(vars)*2)
rownames(tab) <- c(rep(vars, each = 2))
for(i in 1:length(mods)){mod <- mods[[i]]; tabinmod <- unique(rownames(tab)[rownames(tab) %in% rownames(mod)])
  stars <- rep("", length(tabinmod)); stars[mod[tabinmod,4] < .1] <- "^{+}"; stars[mod[tabinmod,4] < .05] <- "^{*}"; 
  stars[mod[tabinmod,4] < .01] <- "^{**}"; stars[mod[tabinmod,4] < .001] <- "^{***}";  
  tab[rownames(tab) %in% rownames(mod),i] <- c(rbind(paste(myround(mod[tabinmod,1],2), stars, sep = ""), paste("(", myround(mod[tabinmod,2],2), ")", sep = "")))}
rownames(tab) <- c(rbind(vars, "")); rownames(tab)[c(1,3,5,7,9,11,13,15,17,19,21,23,25,27)] <- 
  c("Treated (Large firms ben./Small firms harmed)","Age","Male","Black","Latino","AAPI","Other non-white","College-educated","Income","Employed","Unemployed",
    "Party (D=1,R=7)","Ideology (L=1,C=7)","Intercept") 
N <- c(sum(summary(mod1)$df[1:2]),sum(summary(mod2)$df[1:2]),sum(summary(mod3)$df[1:2]),sum(summary(mod4)$df[1:2]))
N <- paste("\\multicolumn{1}{c}{", N, "}", sep = ""); tab <- rbind(tab, N)
addtorow <- list(); addtorow$pos <- list(28,29); addtorow$command <- c(paste0('\\midrule '),paste0('\\midrule '))
ptable(cbind(rownames(tab), tab), "paper/tables/TableB3.tex")

# CI's
LBs <- c(coef(mod1)[2] - qnorm(.975)*summary(mod1)$coef[2,2],
  coef(mod2)[2] - qnorm(.975)*summary(mod2)$coef[2,2],
  coef(mod3)[2] - qnorm(.975)*summary(mod3)$coef[2,2],
  coef(mod4)[2] - qnorm(.975)*summary(mod4)$coef[2,2])
UBs <- c(coef(mod1)[2] + qnorm(.975)*summary(mod1)$coef[2,2],
  coef(mod2)[2] + qnorm(.975)*summary(mod2)$coef[2,2],
  coef(mod3)[2] + qnorm(.975)*summary(mod3)$coef[2,2],
  coef(mod4)[2] + qnorm(.975)*summary(mod4)$coef[2,2])
CIs <- apply(cbind(LBs,UBs), 2, myround, 2)

est <- tab[1,]
ci <- paste("\\multicolumn{1}{c}{", paste("$[", CIs[,1], ",", CIs[,2], "]$\\phantom{'''''}", sep = ""), "}", sep = "")
ci <- gsub("-","\\\\sminus", ci); ci <- gsub("0\\.","\\.", ci);
N <- c(sum(summary(mod1)$df[1:2]),sum(summary(mod2)$df[1:2]),sum(summary(mod3)$df[1:2]),sum(summary(mod4)$df[1:2]))
N <- paste("\\multicolumn{1}{c}{", N, "\\phantom{''''}}", sep = ""); tab <- rbind(tab, N)
tab <- rbind(est, ci, N)
rownames(tab) <- c("Average treatment effect", "ATE 95\\% CI", "N")
addtorow <- list(); addtorow$pos <- list(1); addtorow$command <- c("")
ptable(cbind(rownames(tab), tab), "paper/tables/tableB1_second.tex")

#####################################
## Heterogeneous treatment effects ##
#####################################

## Employmer size
mod1a <- lm(TradeOpenFavorNum ~ I(treatment=="1")*emp_size_big, data = surv); 

## Employmer size
mod1b <- lm(TradeOpenFavorNum ~ I(treatment=="1")*I(emp_size_comb-1), data = surv); 
betas <- mvrnorm(10000, coef(mod1b), vcov(mod1b))
emp_size_seq <- seq(0, 6, by = 1/3)
control_mat <- data.frame(int = 1, treatment = 0, emp_size_comb = emp_size_seq, interaction = 0)
control_sims <- betas %*% t(control_mat)
treated_mat <- data.frame(int = 1, treatment = 1, emp_size_comb = emp_size_seq, interaction = emp_size_seq)
treated_sims <- betas %*% t(treated_mat)

ate.means <- apply(treated_sims - control_sims, 2, mean)
ate.ub <- apply(treated_sims - control_sims, 2, quantile, .975)
ate.lb <- apply(treated_sims - control_sims, 2, quantile, .025)

# plot(ate.means ~ emp_size_seq, type = "l", ylim = c(-.7,.2))
# lines(ate.ub ~ emp_size_seq, type = "l", ylim = c(-.7,.2))
# lines(ate.lb ~ emp_size_seq, type = "l", ylim = c(-.7,.2))

## 
mod2a <- lm(TradeOpenFavorNum ~ I(treatment=="1")*thrm_corp_pos, data = surv); 

## Corporate thermometer
mod2b <- lm(TradeOpenFavorNum ~ I(treatment=="1")*thrm_corp, data = surv); 
betas <- mvrnorm(10000, coef(mod2b), vcov(mod2b))
thrm_corp_seq <- seq(0, 100, by = 5)
control_mat <- data.frame(int = 1, treatment = 0, thrm_corp = thrm_corp_seq, interaction = 0)
control_sims <- betas %*% t(control_mat)
treated_mat <- data.frame(int = 1, treatment = 1, thrm_corp = thrm_corp_seq, interaction = thrm_corp_seq)
treated_sims <- betas %*% t(treated_mat)

ate.means <- apply(treated_sims - control_sims, 2, mean)
ate.ub <- apply(treated_sims - control_sims, 2, quantile, .975)
ate.lb <- apply(treated_sims - control_sims, 2, quantile, .025)

pdf("paper/figures/FigureB1_top.pdf", height = 4, width = 5)

par(mar = c(2.5, 2, 2, 1)) # Set the margin on all sides to 2
plot(ate.means ~ thrm_corp_seq, type = "l", xlim = c(-1,102), ylim = c(-.78,.6), lwd = 2, 
  axes = FALSE, xlab = "", ylab = "")
lines(ate.ub ~ thrm_corp_seq, type = "l", ylim = c(-.72,.4))
lines(ate.lb ~ thrm_corp_seq, type = "l", ylim = c(-.72,.4))
library(scales)

smoother <- density(as.numeric(surv$thrm_corp), kernel = c("gaussian"), width = 2)
lines(smoother$y + (-.76) ~ smoother$x)
jits <- jitter(surv$thrm_corp, factor = 1.5)
segments(x0 = jits, x1 = jits, y0 = -.8, y1 = -.78, col = alpha("black",.2))

mtext(side=2, "Conditional ATE", cex = .6, line = 1, adj = .6)
axis(2, at = c(-.6,-.3,0,.3,.6), labels = c("-.6  ","-.3  "," 0 ",".3 ",".6 "), 
  cex.axis = .5, tck = -.01, padj = 2.75, lwd = .5) 

mtext(side=1, "Views of Corporate America", cex = .6, line = .6, adj = .5)
axis(1, at = c(0,25,50,75,100), labels = c("Negative","","Neutral","","Positive"), 
  cex.axis = .5, tck = -.01, padj = -3.5, lwd = .5) 

mtext(side=3, "Lucid Experiment: Treatment effect conditional on attitude toward corporations", cex = .7, line = .3, adj = .5)

dev.off()

mods <- list(summary(mod2a)$coef, summary(mod2b)$coef, summary(mod1a)$coef, summary(mod1b)$coef)
vars <- c(rownames(mods[[1]])[2:4], rownames(mods[[2]])[3:4], rownames(mods[[3]])[3:4], rownames(mods[[4]])[3:4])
tab <- matrix(data = NA, ncol = length(mods), nrow = length(vars)*2)
rownames(tab) <- c(rep(vars, each = 2))
for(i in 1:length(mods)){mod <- mods[[i]]; tabinmod <- unique(rownames(tab)[rownames(tab) %in% rownames(mod)])
  stars <- rep("", length(tabinmod)); stars[mod[tabinmod,4] < .05] <- "^{*}"; 
  stars[mod[tabinmod,4] < .01] <- "^{**}"; stars[mod[tabinmod,4] < .001] <- "^{***}";  
  tab[rownames(tab) %in% rownames(mod),i] <- c(rbind(paste(myround(mod[tabinmod,1],2), stars, sep = ""), paste("(", myround(mod[tabinmod,2],2), ")", sep = "")))}
rownames(tab) <- c(rbind(vars, "")); rownames(tab)[c(1,3,5,7,9,11,13,15,17)] <- c("Treated (Large firms ben./Small firms harmed)",
    "Positive view of corporations","Treated$\\cdot$ Pos. view of corps.","View of corporations (0-100)","Treated$\\cdot$ View of corps.",
    "Large employer","Treated$\\cdot$ Large employer","Employer size (0-6)","Treated$\\cdot$ Employer size")
N <- c(sum(summary(mod2a)$df[1:2]),sum(summary(mod2b)$df[1:2]),sum(summary(mod1a)$df[1:2]),sum(summary(mod1b)$df[1:2]))
N <- paste("\\multicolumn{1}{c}{", N, "}", sep = ""); tab <- rbind(tab, N)
addtorow <- list(); addtorow$pos <- list(18,19); addtorow$command <- c(paste0('\\midrule '),paste0('\\midrule '))
ptable(cbind(rownames(tab), tab), "paper/tables/TableB6_top.tex")

mod1a <- lm(TradeDealFavorNum ~ I(treatment=="1")*emp_size_big, data = surv); 
mod1b <- lm(TradeDealFavorNum ~ I(treatment=="1")*I(emp_size_comb-1), data = surv); 
mod2a <- lm(TradeDealFavorNum ~ I(treatment=="1")*thrm_corp_pos, data = surv); 
mod2b <- lm(TradeDealFavorNum ~ I(treatment=="1")*thrm_corp, data = surv); 

mods <- list(summary(mod2a)$coef, summary(mod2b)$coef, summary(mod1a)$coef, summary(mod1b)$coef)
vars <- c(rownames(mods[[1]])[2:4], rownames(mods[[2]])[3:4], rownames(mods[[3]])[3:4], rownames(mods[[4]])[3:4])
tab <- matrix(data = NA, ncol = length(mods), nrow = length(vars)*2)
rownames(tab) <- c(rep(vars, each = 2))
for(i in 1:length(mods)){mod <- mods[[i]]; tabinmod <- unique(rownames(tab)[rownames(tab) %in% rownames(mod)])
  stars <- rep("", length(tabinmod)); stars[mod[tabinmod,4] < .05] <- "^{*}"; 
  stars[mod[tabinmod,4] < .01] <- "^{**}"; stars[mod[tabinmod,4] < .001] <- "^{***}";  
  tab[rownames(tab) %in% rownames(mod),i] <- c(rbind(paste(myround(mod[tabinmod,1],2), stars, sep = ""), paste("(", myround(mod[tabinmod,2],2), ")", sep = "")))}
rownames(tab) <- c(rbind(vars, "")); rownames(tab)[c(1,3,5,7,9,11,13,15,17)] <- 
  c("Treated (Large firms ben./Small firms harmed)",
    "Positive view of corporations","Treated$\\cdot$ Pos. view of corps.","View of corporations (0-100)","Treated$\\cdot$ View of corps.",
    "Large employer","Treated$\\cdot$ Large employer","Employer size (0-6)","Treated$\\cdot$ Employer size")
N <- c(sum(summary(mod2a)$df[1:2]),sum(summary(mod2b)$df[1:2]),sum(summary(mod1a)$df[1:2]),sum(summary(mod1b)$df[1:2]))
N <- paste("\\multicolumn{1}{c}{", N, "}", sep = ""); tab <- rbind(tab, N)
addtorow <- list(); addtorow$pos <- list(18,19); addtorow$command <- c(paste0('\\midrule '),paste0('\\midrule '))
ptable(cbind(rownames(tab), tab), "paper/tables/TableB7_top.tex")

###############
## Mediation ##
###############

## Trade openness.
set.seed(12345)
surv$treatment <- as.numeric(surv$treatment == "1")
med.fit1 <- suppressWarnings(polr(as.factor(med_job) ~ treatment + age + gender + race + college + income + employed + unemployed + 
  party + ideology, data = surv, Hess = TRUE))
out.fit1 <- lm(TradeOpenFavorNum ~ med_job + treatment + age + gender + race + college + income + employed + unemployed + 
  party + ideology, data = surv)
med.out1 <- mediate(med.fit1, out.fit1, treat = "treatment", mediator = "med_job",
  robustSE = FALSE, sims = 2000)

med.fit2 <- suppressWarnings(polr(as.factor(med_inq) ~ treatment + age + gender + race + college + income + employed + unemployed + 
  party + ideology, data = surv, Hess = TRUE))
out.fit2 <- lm(TradeOpenFavorNum ~ med_inq + treatment + age + gender + race + college + income + employed + unemployed + 
  party + ideology, data = surv)
med.out2 <- mediate(med.fit2, out.fit2, treat = "treatment", mediator = "med_inq",
  robustSE = FALSE, sims = 2000)

med.fit3 <- suppressWarnings(polr(as.factor(med_pol) ~ treatment + age + gender + race + college + income + employed + unemployed + 
  party + ideology, data = surv, Hess = TRUE))
out.fit3 <- lm(TradeOpenFavorNum ~ med_pol + treatment + age + gender + race + college + income + employed + unemployed + 
  party + ideology, data = surv)
med.out3 <- mediate(med.fit3, out.fit3, treat = "treatment", mediator = "med_pol",
  robustSE = FALSE, sims = 2000)

# 
tabB8topleft <- rbind( 
c(summary(med.out1)$tau.coef, summary(med.out1)$tau.p, summary(med.out1)$tau.ci),

c(summary(med.fit2)$coef[1,c(1,3)], summary(med.fit2)$coef[1,1] - 1.96*summary(med.fit2)$coef[1,2], summary(med.fit2)$coef[1] + 1.96*summary(med.fit2)$coef[1,2]),
c(summary(med.out2)$d1, summary(med.out2)$d1.p, summary(med.out2)$d1.ci),
c(summary(med.out2)$z1, summary(med.out2)$z1.p, summary(med.out2)$z1.ci),

c(summary(med.fit3)$coef[1,c(1,3)], summary(med.fit3)$coef[1,1] - 1.96*summary(med.fit3)$coef[1,2], summary(med.fit3)$coef[1] + 1.96*summary(med.fit3)$coef[1,2]),
c(summary(med.out3)$d1, summary(med.out3)$d1.p, summary(med.out3)$d1.ci),
c(summary(med.out3)$z1, summary(med.out3)$z1.p, summary(med.out3)$z1.ci),

c(summary(med.fit1)$coef[1,c(1,3)], summary(med.fit1)$coef[1,1] - 1.96*summary(med.fit1)$coef[1,2], summary(med.fit1)$coef[1] + 1.96*summary(med.fit1)$coef[1,2]),
c(summary(med.out1)$d1, summary(med.out1)$d1.p, summary(med.out1)$d1.ci),
c(summary(med.out1)$z1, summary(med.out1)$z1.p, summary(med.out1)$z1.ci)
)
tabB8topleft[,c(1,3,4)] <- apply(tabB8topleft[,c(1,3,4)], 2, myround, 2)

sum.med.fit1 <- cbind(summary(med.fit1)$coef, 2*pnorm(-abs(summary(med.fit1)$coef[,3])))
sum.med.fit2 <- cbind(summary(med.fit2)$coef, 2*pnorm(-abs(summary(med.fit2)$coef[,3])))
sum.med.fit3 <- cbind(summary(med.fit3)$coef, 2*pnorm(-abs(summary(med.fit3)$coef[,3])))
mods <- list(sum.med.fit2, summary(out.fit2)$coef, sum.med.fit3, summary(out.fit3)$coef, sum.med.fit1, summary(out.fit1)$coef)
vars <- c(rownames(mods[[1]])[1], rownames(mods[[2]])[2], rownames(mods[[4]])[2], rownames(mods[[6]])[c(2,4:15)])
tab <- matrix(data = NA, ncol = length(mods), nrow = length(vars)*2)
rownames(tab) <- c(rep(vars, each = 2))
for(i in 1:length(mods)){# i <- 1
  mod <- mods[[i]]; tabinmod <- unique(rownames(tab)[rownames(tab) %in% rownames(mod)])
  stars <- rep("", length(tabinmod)); stars[mod[tabinmod,4] < .1] <- "^{+}"; stars[mod[tabinmod,4] < .05] <- "^{*}"; 
  stars[mod[tabinmod,4] < .01] <- "^{**}"; stars[mod[tabinmod,4] < .001] <- "^{***}";  
  tab[rownames(tab) %in% rownames(mod),i] <- c(rbind(paste(myround(mod[tabinmod,1],2), stars, sep = ""), paste("(", myround(mod[tabinmod,2],2), ")", sep = "")))
}
rownames(tab) <- c(rbind(vars, "")); rownames(tab)[c(1,3,5,7,9,11,13,15,17,19,21,23,25,27,29,31)] <- 
  c("Treated (Large firms ben./Small firms harmed)","Mediator: Econ. inequality","Mediator: Socio-pol. inequality","Mediator: Job insecurity",
    "Age","Male","Black","Latino","AAPI","Other non-white","College-educated","Income","Employed","Unemployed",
    "Party (D=1,R=7)","Ideology (L=1,C=7)") 
N <- c(summary(med.fit2)$n,sum(summary(out.fit2)$df[1:2]),summary(med.fit3)$n,sum(summary(out.fit3)$df[1:2]),summary(med.fit1)$n,sum(summary(out.fit1)$df[1:2]))
N <- paste("\\multicolumn{1}{c}{", N, "}", sep = ""); tab <- rbind(tab, N)
addtorow <- list(); addtorow$pos <- list(32,33); addtorow$command <- c(paste0('\\midrule '),paste0('\\midrule '))
ptable(cbind(rownames(tab), tab), "paper/tables/TableB9.tex")

## Trade agreements.
set.seed(12345)
surv$treatment <- as.numeric(surv$treatment == "1")
med.fit1 <- suppressWarnings(polr(as.factor(med_job) ~ treatment + age + gender + race + college + income + employed + unemployed + 
  party + ideology, data = surv, Hess = TRUE))
out.fit1 <- lm(TradeDealFavorNum ~ med_job + treatment + age + gender + race + college + income + employed + unemployed + 
  party + ideology, data = surv)
med.out1 <- mediate(med.fit1, out.fit1, treat = "treatment", mediator = "med_job",
  robustSE = FALSE, sims = 2000)

med.fit2 <- suppressWarnings(polr(as.factor(med_inq) ~ treatment + age + gender + race + college + income + employed + unemployed + 
  party + ideology, data = surv, Hess = TRUE))
out.fit2 <- lm(TradeDealFavorNum ~ med_inq + treatment + age + gender + race + college + income + employed + unemployed + 
  party + ideology, data = surv)
med.out2 <- mediate(med.fit2, out.fit2, treat = "treatment", mediator = "med_inq",
  robustSE = FALSE, sims = 2000)

med.fit3 <- suppressWarnings(polr(as.factor(med_pol) ~ treatment + age + gender + race + college + income + employed + unemployed + 
  party + ideology, data = surv, Hess = TRUE))
out.fit3 <- lm(TradeDealFavorNum ~ med_pol + treatment + age + gender + race + college + income + employed + unemployed + 
  party + ideology, data = surv)
med.out3 <- mediate(med.fit3, out.fit3, treat = "treatment", mediator = "med_pol",
  robustSE = FALSE, sims = 2000)

# 
tabB8topbottom <- rbind( 
c(summary(med.out1)$tau.coef, summary(med.out1)$tau.p, summary(med.out1)$tau.ci),

c(summary(med.fit2)$coef[1,c(1,3)], summary(med.fit2)$coef[1,1] - 1.96*summary(med.fit2)$coef[1,2], summary(med.fit2)$coef[1] + 1.96*summary(med.fit2)$coef[1,2]),
c(summary(med.out2)$d1, summary(med.out2)$d1.p, summary(med.out2)$d1.ci),
c(summary(med.out2)$z1, summary(med.out2)$z1.p, summary(med.out2)$z1.ci),

c(summary(med.fit3)$coef[1,c(1,3)], summary(med.fit3)$coef[1,1] - 1.96*summary(med.fit3)$coef[1,2], summary(med.fit3)$coef[1] + 1.96*summary(med.fit3)$coef[1,2]),
c(summary(med.out3)$d1, summary(med.out3)$d1.p, summary(med.out3)$d1.ci),
c(summary(med.out3)$z1, summary(med.out3)$z1.p, summary(med.out3)$z1.ci),

c(summary(med.fit1)$coef[1,c(1,3)], summary(med.fit1)$coef[1,1] - 1.96*summary(med.fit1)$coef[1,2], summary(med.fit1)$coef[1] + 1.96*summary(med.fit1)$coef[1,2]),
c(summary(med.out1)$d1, summary(med.out1)$d1.p, summary(med.out1)$d1.ci),
c(summary(med.out1)$z1, summary(med.out1)$z1.p, summary(med.out1)$z1.ci)
)
tabB8topbottom[,c(1,3,4)] <- apply(tabB8topbottom[,c(1,3,4)], 2, myround, 2)

sum.med.fit1 <- cbind(summary(med.fit1)$coef, 2*pnorm(-abs(summary(med.fit1)$coef[,3])))
sum.med.fit2 <- cbind(summary(med.fit2)$coef, 2*pnorm(-abs(summary(med.fit2)$coef[,3])))
sum.med.fit3 <- cbind(summary(med.fit3)$coef, 2*pnorm(-abs(summary(med.fit3)$coef[,3])))
mods <- list(sum.med.fit2, summary(out.fit2)$coef, sum.med.fit3, summary(out.fit3)$coef, sum.med.fit1, summary(out.fit1)$coef)
vars <- c(rownames(mods[[1]])[1], rownames(mods[[2]])[2], rownames(mods[[4]])[2], rownames(mods[[6]])[c(2,4:15)])
tab <- matrix(data = NA, ncol = length(mods), nrow = length(vars)*2)
rownames(tab) <- c(rep(vars, each = 2))
for(i in 1:length(mods)){# i <- 1
  mod <- mods[[i]]; tabinmod <- unique(rownames(tab)[rownames(tab) %in% rownames(mod)])
  stars <- rep("", length(tabinmod)); stars[mod[tabinmod,4] < .1] <- "^{+}"; stars[mod[tabinmod,4] < .05] <- "^{*}"; 
  stars[mod[tabinmod,4] < .01] <- "^{**}"; stars[mod[tabinmod,4] < .001] <- "^{***}";  
  tab[rownames(tab) %in% rownames(mod),i] <- c(rbind(paste(myround(mod[tabinmod,1],2), stars, sep = ""), paste("(", myround(mod[tabinmod,2],2), ")", sep = "")))
}
rownames(tab) <- c(rbind(vars, "")); rownames(tab)[c(1,3,5,7,9,11,13,15,17,19,21,23,25,27,29,31)] <- 
  c("Treated (Large firms ben./Small firms harmed)","Mediator: Econ. inequality","Mediator: Socio-pol. inequality","Mediator: Job insecurity",
    "Age","Male","Black","Latino","AAPI","Other non-white","College-educated","Income","Employed","Unemployed",
    "Party (D=1,R=7)","Ideology (L=1,C=7)") 
N <- c(summary(med.fit2)$n,sum(summary(out.fit2)$df[1:2]),summary(med.fit3)$n,sum(summary(out.fit3)$df[1:2]),summary(med.fit1)$n,sum(summary(out.fit1)$df[1:2]))
N <- paste("\\multicolumn{1}{c}{", N, "}", sep = ""); tab <- rbind(tab, N)
addtorow <- list(); addtorow$pos <- list(32,33); addtorow$command <- c(paste0('\\midrule '),paste0('\\midrule '))
ptable(cbind(rownames(tab), tab), "paper/tables/TableB10.tex")


